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ABSTRACT 


Data  projection  (forecasting)  methods  developed  in  the  book  "Time  Series 


Analysis  Forecasting  and  Control"  by  Box  and  Jenkins  (published  by  Holden-Day 
revised  edition  1976)  are  illustrated  using  missile  data  supplied  by  Quality 
Evaluation  Division  of  White  Sands  Missile  Range.  The  prpcess  of  model 
building  making  iterative  use  of  identif ication  (using  the  autocorrelation  func- 
tion) fitting  (maximum  likelihood  estimation)  and  diagnostic  checking  (analysis 
of  residuals)  is  illustrated  in  the  building  of  an  appropriate  stochastic 
difference  equation  model. 


AMS  (MOS)  Subject  Classifications:  60G25,  62M10,  62M20,  62N15 

Key  Words:  Missile,  Trajectory,  Data  projection.  Forecasting  f 
Stochastic  difference  equations 


McmiM  m 


It  is  shown  in  detail  how  all  the  following  may  be  calculated  directly 
from  the  model 

1)  the  projection  (forecast)  function 

2)  the  memory  function 

3)  the  error  of  the  projected  value  at  any  lead  time 

4)  the  updating  of  the  projection  function. 
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ADAPTIVE  MINIMUM  MEAN  SQUARE  ERROR  FORECAST  OF  MISSILE 


TRAJECTORY  USING  STOCHASTIC  DIFFERENCE  EQUATION  MODELS 
G.  E.  P.  Box  and  Lars  Pallesen 

1.  Introduction 

The  objective  of  this  report  is  to  illustrate  data  projection  using 
the  Box  and  Jenkins  text. 

The  data  used  for  illustration  here  were  furnished  by  Mr.  Paul  h.  Thrasher 
of  the  Quality  Evaluation  Division  of  White  Sands  Missile  range. 

2.  Model  Form 


Reasons  are  presented  in  B &J  Chapter  4 for  employing  time  series 
models,  which  are  stochastic  difference  equations  of  the  form 


3.  Identification,  Fitting  and  Checking  of  Model 

The  data  series  we  are  considering  consists  of  246  consecutive 
observations  of  the  x-coordinate  of  a missile  trajectory.  The  observations, 
z^;  t = 1,  2,  . . . , 246,  were  made  with  constant  sampling  interval  and  there 
are  no  missing  or  obviously  aberrant  values. 

Modeling  such  a time  series  is  conceived  of  as  an  iterative  process 
involving  three  stages;  identification,  fitting  and  diagnostic  checking. 
Identification  is  first  performed  along  the  lines  of  Chapter  6 in  B&J  . 

Plotting  the  data  (Figure  la)  shows  a smooth  nonstationary  series, 
whose  autocorrelation  function  (Figure  1(b))  dies  out  extremely  slowly. 

After  differencing  three  times  the  series  appears  stationary  and 

its  sample  autocorrelation  and  partial  autocorrelation  function  (Figures  Ic  and  Id) 
suggest  that  a reasonable  model  for  z^  should  include  a few  moving 
average  parameters  of  low  order.  A clear  identification  is  not  possible  at 
this  point  but  a stochastic  difference  equation  model  of  order  (0,3,3) 

z^  = (1  - OjB  - e^B^  - e^B^a^  (4) 

is  considered  worthy  of  being  tentatively  entertained. 

Fitting  this  model  by  the  method  of  Chapter  7 in  B&J  gives  the 

parameter  estimates,  residual  sum  of  squares  (RSS)  and  the  residual  mean 

square  (RMS)  listed  in  Table  I.  If  this  model  is  adequate  the  RMS  value 

2 2 

provides  an  estimate  of  the  variance,  o-  = E(a  ),  v’hlch  is  the  one 

d t 

step  ahead  forecast  error  variance. 
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Hgure  lb.  Sample  aufocorrclallon  function  of  original  scries 
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Figure  Ic.  Sample  autocorrelation  function  for 
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(Pfd,q) 


,2,2 


= (1  - 0^6-028^)3^ 


0^  = .716 


'^Moduli  of  roots 


ots:  1.  39;  1.  39\ 
i.  e.  stable  j 


V^z^  = (1  - 0jB-02B^-02B^)a^ 


0^  = .662 


02  =-''4 


03  = -.425 


/Moduli  of  roots 


l:ll 

(::  54 

)ts;  1.  14;  1.  14;  1.  83\ 
i.  e.  stable  ) 


V z^  = (1-0^B- 

■®2®  -636  )a^ 

0^  = 1.  731 

ri.75 

W.72 

0-  = -.  776 

r-.76 

2 

1 -.79 

0,  = -. 104 

f -.08 

3 

t^-.13 

^Moduli  of  roots: 

1.  014;  1.  OH;  9.  39’ 

, i.  e. 

stable 

RMS  (DF 


1.03  (240) 


Table  I Continued 


(P,d,q) 

Mjdel 

RSS 

RMS 

(DF) 

(0,  3,4) 

= (l-F,B-e-B^-0,B^-e  B‘*)a^ 
t 1 2 3 4 t 

203. 

. 85 

(239) 

0^  = 1.938 

1. 99 

1. 89 

©2  = -1. 030  1 

r 

- .95 

-1. 11 

03  = - . 146  1 

f - . 02 

1 - .27 

©4  = . 173  < 

r 

. 09 

(Moduli  of  roots : 

; 1.  13;  1.  13; 1.  59; 

2.  86) 

(0,3,5) 

= (l-0,B-0^B^-0  B^-0  b'^-0^B^) 
t i Z 3 4 D 

192. 

. 81 

(238) 

Sj  = 2.078  1 

r 2. 11 

2.  04 

% = -'•29'  i 

r -1.21 

1 -1. 37 

e,  .-,.115  1 

f . 00 

1 - . 23 

. . 395  1 

[ .51 

. 28 

0^  = - . 131  < 

0 

r - . 04 
[ - . 22 

(Moduli  of  roots: 

1.11;1. 11;1.79;1.  79 

1.  93) 

(0,  3,  6) 

V^z^  = (1-0,B-0^B,-0,B^-0  b‘*-0^B® 
t 1 Z 3 4 D 

192. 

. 81 

(237) 

(roots  o,  k. ) 

0,  a;  0 

Diagnostic  checking  (Chapter  8 in  B&J)  involves  examination  of 
the  residuals  (the  estimated  a^'s)  left  after  fitting  this  model  to  seek  for 
departures  from  the  "white  noise"  form.  One  way  of  doing  this  is  to  sub- 
mit the  residual  a^  sequence  to  the  identification  procedure  previously 
3 

applied  to  V . In  fact  the  autocorrelation  function  of  the  residuals 
a^'s  , Figure  2(a),  suggests  that  while  most  of  the  dependence  is  being 
accounted  for  by  the  model,  some  significant  low  order  autocorrelation 
remains,  indicating  some  additional  0 parameters  are  needed.  Notice, 
that  the  diagnostic  checking  of  the  model  (4)  reveals  model  inadequacy  and 
also  identifies  in  which  way  the  model  should  be  modified. 

After  another  cycle  the  (0,3,!5)  moae. 

= (1  - e^B  - e^B^  - e^B^  - b^b"^  - e^B^)a^  (5) 

is  considered,  and  it  fits  the  data  very  well,  leaving  residuals,  Figure  2(b), 
which  look  like  white  noise.  Figure  2(c)  shows  the  sample  autocorrelations 
of  the  residuals.  This  fitted  model  along  with  some  other  contenders  are 
listed  in  Table  I.  Additional  models  are  fitted  as  a check  that  additional 
parameters  would  not  substantially  improve  matters  (overfitting),  and  also 
to  demonstrate  that  the  chosen  number  of  differencings  is  appropriate. 

4.  Checking  9(B) 

Regarding  the  operator 

0(B)  = 1 - 0B  - - 0^B^  - B^b'^  - B^B^  (6) 

as  a polynomial  in  B,  Lt  is  shown  in  B &J  that  a necessary  requirement 
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Figure  2a.  Autocorrelation  function  of  residuals  from  the  (0,  3,  3)  model 
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Figure  2c.  Autocorrelation  function  of  residuals  from  the  (0,  3,  5)  model 


-12- 


for  a sensible  model  is  that  the  zeroes  of  this  polynomial  Pe  outside  the 
unit  circle  (invertibility  property). 

It  is  important  to  check  this  and  the  moduli  of  the  roots  given 
in  Table  I indicate  that  the  model  is  indeed  invertible. 


5.  Forecasts 

Accepting  that  the  (0,  3,  5)  model  provides  an  adequate  representa- 
tion of  the  system  (with  the  (0,  3,4)  model  as  a close  runner-up)  the  forecasts 
produced  are  most  easily  calculated  from  the  difference  equation  itself 
(see  Chapter  5 of  B&J).  From  Equation  (5)  we  find 


"t  = '"t-1  ' '"t-2  ^ "t-3 


+ a^  - eia^_j  - ' 63^.3  ' 64^-4  " ®5^-5  ’ 


(7) 


Then  by  taking  conditional  expectations  of  ^t+2’ ' ’ ' ’ ^t+f 

t (as  described  in  B&I  p.  130)  the  1,2,3,...,!,...  step  ahead  forecasts 


are; 


z^(l)  = 3z^  - 3z^_j  + z^_2  - - e^a^.i  - e3a^_2  - e^a^.3  ' 05^.4 

z^(Z)=  3z^(l)  - 3z^  + - e2a^  - 03aj_j  - e^a^,^  ' Vt-3 

S^(3)  = 3z^(2)  - 3z^(l)  + - 633^  - e^a^.i  - 659^2 

S^(4)  = 3z^(3)  - 3z^(2)  +z(l)  - e^a^  - e5a^_^ 

2^(5)  = 3z^4)  - 3z^(3)  - z(2)  - O^a^ 
z (f)  = 3z  (i-1)  - 3z(l-2)  - z(l-3) 

• L ^ 


(8) 


I > 6 


In  practice  of  course  this  is  done  automatically  by  the  computer. 
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Table  II 


Obs  # 

Actual  value 

Model 
(0,  3,  5) 

Forecasts 
Model 
(0,  3,4) 

Model 
(0,  3,  3) 

201 

13225. 08 

13224. 78 

13224. 80 

13224. 99 

202 

13306.74 

13305. 80 

13305.  94 

13306.46 

203 

13387.51 

13386.  70 

13386. 77 

13387. 51 

204 

13468.42 

13467.  20 

13467.  23 

13468.  14 

205 

13549. 74 

13547. 33 

13547. 34 

'T'3548.  34 

206 

13628. 61 

13627.  10 

13627.  08 

13628.  12 

207 

1 

13708. 78 

13706.49 

13706.46 

13707.  48 

208 

13788. 67 

13785. 52 

13785.48 

13786.40 

209 

13868. 21 

13864. 18 

13864.14 

13864. 91 

210 

13947. 30 

13942.47 

13942.44 

13942. 99 

14- 


For  illustration,  the  forecasts  produced  by  this  model  with  an  origin 
(for  all  forecasts)  at  t = ZOO  is  shown  in  Figure  3.  It  will  be  noticed 


that  the  forecasts  are  in  very  close  agreement  with  the  actual  values. 

Even  the  10  step  ahead  forecast  is  hardly  distinguishable  from  the  actually 
observed  value. 

Table  II  lists  the  actual  values  and  the  forecasts  numerically.  The 
forecasts  produced  by  the  models  (0,  3,4)  and  (0,  3,  3)  are  also  very  good 
and  they  are  included  for  comparison. 

6.  Error  of  Forecasts 

In  order  to  determine  the  error  of  the  forecasts,  it  is  helpful  to  write 
the  model  (1)  in  random  shock  form.  Thus  formally 

e (B) 

z.  = -2 a.  = 4j(B)a  (9) 

V 4.^(8) 

where 

2 

4>(B)  = 1 + + 4^2®  + • • • (10) 

And  it  is  shown  in  B&J  p.  126-128  that  the  lead  f forecast  error  is 


L6 


(13) 


• " ‘•'/-I't  + l 


V^ence  the  variance  of  the  forecast  error  is 


var[e^(f )]  = 

2 2 2 2 

= (1  + 4jj  + 4;^  + • • • ■•■  .l^  *^3  ■ 

For  the  fitted  model  (5)  the  4j-weights  are  calculated  by  equating  coefficients 
in  (15),  B &J  pp.  132-134. 

(1  - 3B  + 3B^  - B^)(l  + 4)jB  + 4^2^^  + ...)  = (1  - BjB  - - Q^B^  - 0^B‘*  - 0^B^) 

(15) 

Specifically  we  find  the  4.^  values  given  in  Table  III.  Using  the  esti- 
mated = . 81  from  Table  I,  the  variance  of  the  forecast  error  is  given 

for  I = 1,  2,  . . . , 10  . The  last  column  in  Table  III  lists  ± 2 standard 
errors,  corresponding  to  approximately  95%  probability  intervals  for  the 
forecasts.  We  note,  that  these  probability  intervals  are  so  narrow,  that 
they  cannot  be  distinguished  from  the  forecasts  themselves  in  a plot  like 
Figure  3. 

All  that  is  needed  to  compute  forecasts  and  the  standard  deviations 
of  forecast  errors,  is  given  here.  What  appears  in  the  following  Appendices 
is  not  necessary  for  calculation,  but  does  illuminate  the  nature  of  the 
projection  process. 
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Appendix  A - Integral  forms 


As  discussed  in  Chapter  4 pp.  103-114  of  B &cj,  the  equivalent 
integrated  forn.  of  the  model  of  Equation  (5),  is  of  some  interest  also. 

In  this  form  the  observations  appear  as  a linear  aggregates  of  past 
random  shocks,  their  difference,  sum,  sum  of  sums,  etc,  , plus  a new 
random  shock.  Specifically  the  Integrated  model  form 

degenerates  to  different  models  from  Table  I when  certain  of  the  \-co- 
efficients  are  taken  to  be  zero.  Table  IV  links  models  from  Table  I to 
their  equivalent  Integrated  forms,  and  lists  estimated  \ coefficients 
which  can  be  calculated  from  the  estimated  B's  . Conversion  formulas 
for  the  models  under  consideration  are  given  in  Table  V,  but  can  more 
generally  be  found  from  equating  coefficients  in  Equation  4.  3.  21,  p.  112 
in  B&J. 
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Table  IV 


integrated  model  forms 


Model 

^2 

^l 

^0 

s 

^2 

RSS 

(0,2,2) 

- 

- 

. 483 

. 800 

- 

410. 

(0,2,  3) 

- 

.425 

. 035 

. 878 

- 

348. 

(0,  3,  3) 

- 

- 

1.  104 

. 016 

. 149 

247. 

(0,  3,4) 

. 173 

. 627 

. 197 

. 06  5 

203. 

(0,  3,  5) 

. 131 

129 

.716 

. 140 

. 064 

192. 

Table  V 


Conversion  formulae  6 to  \ 


Model  Formulae 


(0,2,2) 

Kq  = 1 . 

(0,2,  3) 

^1  = -®3 

"O  = 1 + 02  + 203 

N =1-01-02-03 

■ 

(0,  3,  3) 

S = ’ - ®3 

Xj  = 1 + + 263 

^2  =1-01-0,-63 

(0,  3,4) 

^ = 1 - O3  - 36, 

= 1 + 0^  + 20  + 30, 

1 2 3 4 

^2  =1-01-0,-63-6, 

(0,  3,5) 

^2=-^ 

X , = 0.  + 40 

X^  = 1 - 63  - 30,  - 60^ 

X,  = 1 + 0,  + 20  + 30^  + 40 

1 2 3 4 5 

Xz  = 1 - 0,  - O2  - ^ - O5  - ®5 

Appendix  B - The  eventual  forecast  function 


One  question  of  interest  is  what  function  is  being  selected  for  pro- 
jecting the  forecasts,  i.  e.  what  is  the  forecast  function.  It  is  shown  in 
B&J  p.  139  that  depending  on  the  nature  of  the  left  hand  operator,  the 
model  (1)  could  call  for  forecasts  lying  on  an  updating  function  that  could 
consist  of  any  combination  of  polynomials,  exponentials  and  sine  and 
cosine  waves.  What  forecast  function  does  the  model  imply  for  the 
present  fitted  (0,  3,  5)  model? 

The  eventual  forecast  function  for  the  (0,  3,  5)  model  satisfies  the 
difference  equation 


V z^(jf)  = 0 


(B-1) 


nd 


which  has  as  its  solutions  a polynomial  in  f of  2 degree 


2,(0  = b|,'' . b|*>, 


(B-^) 


and  applies  for  i > q - p - d (i.  e.  I > 2)  . 

In  other  words  the  model  (0,  3,  5)  implies,  that  the  forecasted 
future  values  from  any  time  origin  t will,  except  for  slight  deviations 
at  the  first  two  lead-times,  follow  a quadratic  curve.  (The  (0,  3,4) 
model  which  fits  slightly  less  well  Implies  that  only  one  initial  deviation 
occurs,  while  the  (0,  3,  3)  model  implies  that  all  forecasts  lie  on  a 
quadratic  curve). 

Although  the  forecasts  are  best  calculated  directly  from  the  difference 


equation  as  above  it  is  enlightening  to  further  consider  their  nature. 


As  the  origin  of  forecasts  is  advanced  the  calculating  process  re- 
quires that  coefficients  b^,  and  b^  are  sequentially  updated.  For 
example  the  updating  formulae  for  the  (0,  3,  5)  model  can  be  found 
directly  by  relating  (B-2)  to  the  forecasting  formula  from  the  integrated 
model. 

We  find  that  the  updating  formulae  derived  below  are 


Note  that  the  first  terms  on  the  right  of  (B-3)  simply  allow  for  movement 
of  the  origin  without  changing  the  polynomial.  The  term  involving  the 
last  random  shock  a^  appropriately  updates  the  coefficient. 

The  updating  formulae  (B-3)  are  derived  as  follows.  We  have 
from  Equation  (A-1)  in  Appendix  A that 


= X.  _Va^  .,fX.,a  , + \„Sa  , 

t+f  -2  t+i-1  -1  t+i-1  0 t+jC-1 


+ \,  S a - + X._  S a , + a 
1 t+/-l  2 t+i-1  t+i 


(B-4) 


Assuming  f > 2 and  taking  expectations  at  origin  t we  find 
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N > 


2 3 

ri)  = E(>.  Sa,  . ,)+E(\,S  a ,)  + E(\^S  a ,) 

' ' 0 t+f -r  ' 1 t+i -r  ' 2 t+f-1 


~ ( ^0  ^ S a^  j ^ Xj  S a^) 

+ ^2  5=,) 

= (Xq  Sa^  + X^  S a^_^  ^2  ® ®t-l  ^ ^2  ^ ^t-2^ 

+ ^(X^Sa^  + + j X^  S a^) 


+ t (]  X^  S a^) 


(B-5) 


The  coefficients  b^,  , b^  in  Equation  (B-2)  are  now  identified  as 


'’o'  = S ® “t  ^ ’'l  ®t-l  * ^2  °t-l*  ^2^^  “t-2 


bf  = X,Sa,  tX2s"aj_j  + iX,Sat 


U^(t)  _ 1 \ c = 
^2  - ^ ^2  ^ ®t 


(B-6) 


Now  it  is  seen  that  (B-6)  can  be  rewritten  as  (B-3) 
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Appendix  C - How  are  the  cata  used  in  the  forecast? 


Still  another  way  to  interpret  the  forecasts  is  as  a weighted  sum  of 
previous  observations;  Writing  (5)  as 

— Zj  = T,,B)Zj  = a^  (C-1) 

where 

■n-(B)  = 1 - tTjB  - . (C-2) 

we  find  that 

z^(t)  = Tr^z^(i-l)  + + ...  (C-3) 

where  z^(-h)  is  taken  to  mean  z^  ^ for  h = 0, 1,  2,  . . . . The  t7 -weights 
can  be  found  by  equating  coefficients  in  the  following  identity  after  the 
0-estimates  have  been  substituted* 

= e(B)  tt(B) 

(1  - 3B  + 3B^  - B^)  = (1  - GjB-e^B^  - G^B^  - G^B^  - G^B^Kl  - tt^B  - rr^B^  - tt^B^  - . . . ) 

(C-4) 

The  iT-weights  (also  denoted  by  are  given  in  Figure  4;  thus 

for  example 

z^(l)  = . 922  X + . 207  X z^  ^ + . 355  X z^  ^ ' • 039  X z^  ^ 

-.068  Xz^  -,171xz  -.149XZ  , -.143  Xz 

t-4  t-5  t-o  t-7 

-.  107  XZ^_g  - . 079  Xz^_^  - . 046  XZ^_^q  - . 018  Xz^_^^  (^-5) 

+ .008Xz^  _+.027  Xz^  ,_+.04lXz  +.048  XZ 

t-lZ  1-13  t-14  t-lb 

* ■ “'*«  X^.16  * ■ *■■■ 
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-weights  for  th 


The  two  step  ahead  forecast  can  be  found  similarly  by  replacing  by 

z^(l)  and  z^  ^ by  z^  , ,,  and  so  on  for  forecasts  with  higher  lead 
t t-j  t-]+r 

times.  However  these  forecasts  may  also  be  expressed  directly  as 
weighted  sums  of  the  observations  z^,  1’  ^t  2’  ' The  weights 

and  corresponding  to  the  two  and  three  step  ahead  forecasts 

respectively,  are  also  show  in  Figure  4.  In  general  these  weights  may  be 
found  from  the  u^^^-weight  by  means  of  the  formula  (p.  142  Equation  5.  3.  9) 


] 


I -1 


y IT,  it; 
j+i-1  h j 


(jf-h) 


(C-6) 
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